home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / linalg / symmtd.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  6.6 KB  |  229 lines
  1. /* linalg/sytd.c
  2.  * 
  3.  * Copyright (C) 2001 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Factorise a symmetric matrix A into
  21.  *
  22.  * A = Q T Q'
  23.  *
  24.  * where Q is orthogonal and T is symmetric tridiagonal.  Only the
  25.  * diagonal and lower triangular part of A is referenced and modified.
  26.  *
  27.  * On exit, T is stored in the diagonal and first subdiagonal of
  28.  * A. Since T is symmetric the upper diagonal is not stored.
  29.  *
  30.  * Q is stored as a packed set of Householder transformations in the
  31.  * lower triangular part of the input matrix below the first subdiagonal.
  32.  *
  33.  *  * The full matrix for Q can be obtained as the product
  34.  *
  35.  *       Q = Q_N .. Q_2 Q_1
  36.  *
  37.  * where 
  38.  *
  39.  *       Q_i = (I - tau_i * v_i * v_i')
  40.  *
  41.  * and where v_i is a Householder vector
  42.  *
  43.  *       v_i = [1, A(i+2,i), A(i+3,i), ... , A(N,i)]
  44.  *
  45.  * This storage scheme is the same as in LAPACK.  See LAPACK's
  46.  * ssytd2.f for details.
  47.  *
  48.  * See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 8.3 */
  49.  
  50. #include <config.h>
  51. #include <stdlib.h>
  52. #include <gsl/gsl_math.h>
  53. #include <gsl/gsl_vector.h>
  54. #include <gsl/gsl_matrix.h>
  55. #include <gsl/gsl_blas.h>
  56.  
  57. #include <gsl/gsl_linalg.h>
  58.  
  59. int 
  60. gsl_linalg_symmtd_decomp (gsl_matrix * A, gsl_vector * tau)  
  61. {
  62.   if (A->size1 != A->size2)
  63.     {
  64.       GSL_ERROR ("symmetric tridiagonal decomposition requires square matrix",
  65.                  GSL_ENOTSQR);
  66.     }
  67.   else if (tau->size + 1 != A->size1)
  68.     {
  69.       GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
  70.     }
  71.   else
  72.     {
  73.       const size_t N = A->size1;
  74.       size_t i;
  75.   
  76.       for (i = 0 ; i < N - 2; i++)
  77.         {
  78.           gsl_vector_view c = gsl_matrix_column (A, i);
  79.           gsl_vector_view v = gsl_vector_subvector (&c.vector, i + 1, N - (i + 1));
  80.           double tau_i = gsl_linalg_householder_transform (&v.vector);
  81.           
  82.           /* Apply the transformation H^T A H to the remaining columns */
  83.  
  84.           if (tau_i != 0.0) 
  85.             {
  86.               gsl_matrix_view m = gsl_matrix_submatrix (A, i + 1, i + 1, 
  87.                                                         N - (i+1), N - (i+1));
  88.               double ei = gsl_vector_get(&v.vector, 0);
  89.               gsl_vector_view x = gsl_vector_subvector (tau, i, N-(i+1));
  90.               gsl_vector_set (&v.vector, 0, 1.0);
  91.               
  92.               /* x = tau * A * v */
  93.               gsl_blas_dsymv (CblasLower, tau_i, &m.matrix, &v.vector, 0.0, &x.vector);
  94.  
  95.               /* w = x - (1/2) tau * (x' * v) * v  */
  96.               {
  97.                 double xv, alpha;
  98.                 gsl_blas_ddot(&x.vector, &v.vector, &xv);
  99.                 alpha = - (tau_i / 2.0) * xv;
  100.                 gsl_blas_daxpy(alpha, &v.vector, &x.vector);
  101.               }
  102.               
  103.               /* apply the transformation A = A - v w' - w v' */
  104.               gsl_blas_dsyr2(CblasLower, -1.0, &v.vector, &x.vector, &m.matrix);
  105.  
  106.               gsl_vector_set (&v.vector, 0, ei);
  107.             }
  108.           
  109.           gsl_vector_set (tau, i, tau_i);
  110.         }
  111.       
  112.       return GSL_SUCCESS;
  113.     }
  114. }  
  115.  
  116.  
  117. /*  Form the orthogonal matrix Q from the packed QR matrix */
  118.  
  119. int
  120. gsl_linalg_symmtd_unpack (const gsl_matrix * A, 
  121.                           const gsl_vector * tau,
  122.                           gsl_matrix * Q, 
  123.                           gsl_vector * diag, 
  124.                           gsl_vector * sdiag)
  125. {
  126.   if (A->size1 !=  A->size2)
  127.     {
  128.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  129.     }
  130.   else if (tau->size + 1 != A->size1)
  131.     {
  132.       GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
  133.     }
  134.   else if (Q->size1 != A->size1 || Q->size2 != A->size1)
  135.     {
  136.       GSL_ERROR ("size of Q must match size of A", GSL_EBADLEN);
  137.     }
  138.   else if (diag->size != A->size1)
  139.     {
  140.       GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
  141.     }
  142.   else if (sdiag->size + 1 != A->size1)
  143.     {
  144.       GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
  145.     }
  146.   else
  147.     {
  148.       const size_t N = A->size1;
  149.  
  150.       size_t i;
  151.  
  152.       /* Initialize Q to the identity */
  153.  
  154.       gsl_matrix_set_identity (Q);
  155.  
  156.       for (i = N - 2; i > 0 && i--;)
  157.     {
  158.           gsl_vector_const_view c = gsl_matrix_const_column (A, i);
  159.           gsl_vector_const_view h = gsl_vector_const_subvector (&c.vector, i + 1, N - (i+1));
  160.           double ti = gsl_vector_get (tau, i);
  161.  
  162.           gsl_matrix_view m = gsl_matrix_submatrix (Q, i + 1, i + 1, N-(i+1), N-(i+1));
  163.  
  164.           gsl_linalg_householder_hm (ti, &h.vector, &m.matrix);
  165.         }
  166.  
  167.       /* Copy diagonal into diag */
  168.  
  169.       for (i = 0; i < N; i++)
  170.         {
  171.           double Aii = gsl_matrix_get (A, i, i);
  172.           gsl_vector_set (diag, i, Aii);
  173.         }
  174.  
  175.       /* Copy subdiagonal into sd */
  176.  
  177.       for (i = 0; i < N - 1; i++)
  178.         {
  179.           double Aji = gsl_matrix_get (A, i+1, i);
  180.           gsl_vector_set (sdiag, i, Aji);
  181.         }
  182.  
  183.       return GSL_SUCCESS;
  184.     }
  185. }
  186.  
  187. int
  188. gsl_linalg_symmtd_unpack_T (const gsl_matrix * A, 
  189.                             gsl_vector * diag, 
  190.                             gsl_vector * sdiag)
  191. {
  192.   if (A->size1 !=  A->size2)
  193.     {
  194.       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
  195.     }
  196.   else if (diag->size != A->size1)
  197.     {
  198.       GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
  199.     }
  200.   else if (sdiag->size + 1 != A->size1)
  201.     {
  202.       GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
  203.     }
  204.   else
  205.     {
  206.       const size_t N = A->size1;
  207.  
  208.       size_t i;
  209.  
  210.       /* Copy diagonal into diag */
  211.  
  212.       for (i = 0; i < N; i++)
  213.         {
  214.           double Aii = gsl_matrix_get (A, i, i);
  215.           gsl_vector_set (diag, i, Aii);
  216.         }
  217.  
  218.       /* Copy subdiagonal into sdiag */
  219.  
  220.       for (i = 0; i < N - 1; i++)
  221.         {
  222.           double Aij = gsl_matrix_get (A, i+1, i);
  223.           gsl_vector_set (sdiag, i, Aij);
  224.         }
  225.  
  226.       return GSL_SUCCESS;
  227.     }
  228. }
  229.